/*
 * Copyright 1990-2008, Mark Little, University of Newcastle upon Tyne
 * and others contributors as indicated 
 * by the @authors tag. All rights reserved. 
 * See the copyright.txt in the distribution for a
 * full listing of individual contributors. 
 * This copyrighted material is made available to anyone wishing to use,
 * modify, copy, or redistribute it subject to the terms and conditions
 * of the GNU Lesser General Public License, v. 2.1.
 * This program is distributed in the hope that it will be useful, but WITHOUT A 
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 
 * PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
 * You should have received a copy of the GNU Lesser General Public License,
 * v.2.1 along with this distribution; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 
 * MA  02110-1301, USA.
 * 
 * (C) 1990-2008,
 */
 /*
 * Copyright (C) 1994-1998,
 *
 * Department of Computing Science,
 * The University,
 * Newcastle upon Tyne,
 * UK.
 *
 * $Id: Random.n,v 1.3 1998/08/28 14:19:53 nmcl Exp $
 */

#include <math.h>

#if defined(NO_INLINES) && !defined(RANDOM_CC_)

#else

#ifndef NO_INLINES
#  define INLINEF inline
#else
#  define INLINEF
#endif

INLINEF double RandomStream::MGen ()
{
    // A multiplicative generator, courtesy I. Mitrani 1992,
    // private correspondence
    // Y[i+1] = Y[i] * 5^5 mod 2^26
    // period is 2^24, initial seed must be odd

    const long int two2the26th = 67108864;	// 2**26

    MSeed = (MSeed * 25) % two2the26th;
    MSeed = (MSeed * 25) % two2the26th;
    MSeed = (MSeed * 5) % two2the26th;

    return (double) MSeed / (double) two2the26th;
}

INLINEF double RandomStream::Uniform () 
{
    // A linear congruential generator based on the algorithm from
    // "Algorithms", R. Sedgewick, Addison-Wesley, Reading MA, 1983.
    // pp. 36-38.

    const long m=100000000;
    const long b=31415821;
    const long m1=10000;

    // Do the multiplication in pieces to avoid overflow
    
    long    p0 = LSeed%m1,
	    p1 = LSeed/m1,
	    q0 = b%m1,
	    q1 = b/m1;

    LSeed = (((((p0*q1+p1*q0)%m1)*m1+p0*q0)%m) + 1) % m;

    // The results of the LC generator are shuffled with
    // the multiplicative generator as suggested by
    // Maclaren and Marsaglia (See Knuth Vol2, Seminumerical Algorithms)

    long choose = LSeed % (sizeof(series)/sizeof(double));
    double result = series[choose];
    
    series[choose] =  MGen();

    return result;
}

INLINEF double RandomStream::Error ()
{
    const long r=100;
    const long N=100*r;
    long i, f[r];
    for (i=0; i<r; i++) f[i]=0;
    for (i=0; i<N; i++) f[(int) (Uniform()*r)]++;
    long t=0;
    for (i=0; i<r; i++) t += f[i]*f[i];
    double rt = (double) r*t;
    double rtN = rt / (double) N - (double) N;
    return 1.0 - (rtN / r);
}

INLINEF double UniformStream::operator() ()
{
    return lo+(range*Uniform());
}

INLINEF Boolean Draw::operator() ()
{
    return (Boolean) (s() >= prob);
}

INLINEF double ExponentialStream::operator() ()
{
    return -Mean*log(Uniform());
}

INLINEF double ErlangStream::operator() ()
{
    double z=1.0;
    for (int i=0; i<k; i++) z*=Uniform();
    return -(Mean/k)*log(z);
}

INLINEF double HyperExponentialStream::operator() ()
{
    double z = (Uniform()>p) ? Mean/(1.0-p) : Mean/p;
    return -0.5*z*log(Uniform());
}

INLINEF double NormalStream::operator() ()
{
    // Use the polar method, due to Box, Muller and Marsaglia
    // Taken from Seminumerical Algorithms, Knuth, Addison-Wesley, p.117

    double X2;

    if (z!=0.0) {
	X2 = z;
	z = 0.0;
    } else {
	double S, v1, v2;
	do {
	    v1 = 2.0*Uniform()-1.0;
	    v2 = 2.0*Uniform()-1.0;
	    S = v1*v1 + v2*v2;
	} while (S>=1.0);

	S = sqrt((-2.0*log(S))/S);
	X2 = v1*S;
	z  = v2*S;
    }

    return Mean + X2*StandardDeviation;
}

#ifdef INLINEF
#  undef INLINEF
#endif

#endif
